Collaborators: Shane Blowes, Jon Chase, Helmut Hillebrand, Michael Burrows, Amanda Bates, Uli Brose, Benoit Gauzens, Laura Antao Assistance: Katherine Lew, Josef Hauser
library(data.table) # for handling large datasets
library(ggplot2) # for some plotting
#library(lme4)
library(nlme) # for ME models
library(beanplot) # for beanplots
library(maps) # for map
library(ggeffects) # marginal effect plots
library(gridExtra) # to combine ggplots together
library(grid) # to combine ggplots together
library(gridExtra)
options(width=500) # turn off most text wrapping
# tell RStudio to use project root directory as the root for this notebook. Needed since we are storing code in a separate directory.
knitr::opts_knit$set(root.dir = rprojroot::find_rstudio_root_file())
# Turnover and covariates assembled by turnover_vs_temperature_prep.Rmd
trends <- fread('output/turnover_w_covariates.csv.gz')
# set realm order
trends[, REALM := factor(REALM, levels = c('Freshwater', 'Marine', 'Terrestrial'), ordered = FALSE)]
# group Marine invertebrates/plants in with All
trends[, taxa_mod2 := taxa_mod]
trends[taxa_mod == 'Marine invertebrates/plants', taxa_mod2 := 'All']
Log-transform some variables, then center and scale.
trends[, tempave.sc := scale(tempave)]
trends[, tempave_metab.sc := scale(tempave_metab)]
trends[, seas.sc := scale(seas)]
trends[, microclim.sc := scale(log(microclim))]
trends[, temptrend.sc := scale(temptrend)]
trends[, temptrend_abs.sc := scale(log(abs(temptrend)))]
trends[, npp.sc := scale(log(npp))]
trends[, mass.sc := scale(log(mass_mean_weight))]
trends[, speed.sc := scale(log(speed_mean_weight+1))]
trends[, lifespan.sc := scale(log(lifespan_mean_weight))]
trends[, thermal_bias.sc := scale(thermal_bias)]
trends[, consumerfrac.sc := scale(consfrac)]
trends[, endothermfrac.sc := scale(endofrac)]
trends[, nspp.sc := scale(log(Nspp))]
trends[, human.sc := scale(human)]
# histograms to examine
cexmain = 0.6
par(mfrow = c(3,5))
invisible(trends[, hist(tempave.sc, main = 'Environmental temperature (°C)', cex.main = cexmain)])
invisible(trends[, hist(tempave_metab.sc, main = 'Metabolic temperature (°C)', cex.main = cexmain)])
invisible(trends[, hist(seas.sc, main = 'Seasonality (°C)', cex.main = cexmain)])
invisible(trends[, hist(microclim.sc, main = 'log Microclimates (°C)', cex.main = cexmain)])
invisible(trends[, hist(temptrend.sc, main = 'Temperature trend (°C/yr)', cex.main = cexmain)])
invisible(trends[, hist(temptrend_abs.sc, main = 'log abs(Temperature trend) (°C/yr)', cex.main = cexmain)])
invisible(trends[, hist(mass.sc, main = 'log Mass (g)', cex.main = cexmain)])
invisible(trends[, hist(speed.sc, main = 'log Speed (km/hr)', cex.main = cexmain)])
invisible(trends[, hist(lifespan.sc, main = 'log Lifespan (yr)', cex.main = cexmain)])
invisible(trends[, hist(consumerfrac.sc, main = 'Consumers (fraction)', cex.main = cexmain)])
invisible(trends[, hist(endothermfrac.sc, main = 'Endotherms (fraction)', cex.main = cexmain)])
invisible(trends[, hist(nspp.sc, main = 'log Species richness', cex.main = cexmain)])
invisible(trends[, hist(thermal_bias.sc, main = 'Thermal bias (°C)', cex.main = cexmain)])
invisible(trends[, hist(npp.sc, main = 'log Net primary productivity', cex.main = cexmain)])
invisible(trends[, hist(human.sc, main = 'Human impact score', cex.main = cexmain)])
panel.cor <- function(x, y, digits = 2, prefix = "", cex.cor, ...)
{
usr <- par("usr"); on.exit(par(usr))
par(usr = c(0, 1, 0, 1))
r <- cor(x, y, use = 'pairwise.complete.obs')
txt <- format(c(r, 0.123456789), digits = digits)[1]
txt <- paste0(prefix, txt)
if(missing(cex.cor)) cex.cor <- 0.8/strwidth(txt)
text(0.5, 0.5, txt) #, cex = cex.cor * r)
}
pairs(formula = ~ REALM + tempave.sc + tempave_metab.sc + seas.sc + microclim.sc + temptrend.sc + temptrend_abs.sc + mass.sc + speed.sc + lifespan.sc + consumerfrac.sc + endothermfrac.sc + nspp.sc + thermal_bias.sc + npp.sc + human.sc, data = trends, gap = 1/10, cex = 0.2, col = '#00000022', lower.panel = panel.cor)
Mass and lifespan look tightly correlated, but r only 0.56…? Tempave_metab and lifespan don’t look tightly correlated, but r= -0.81 Tempave_metab and speed don’t look tightly correlated, but r= -0.83 Lifespan and speed don’t look tightly correlated, but r = 0.73
Just turnover
cat('Overall # time-series: ', nrow(trends), '\n')
Overall # time-series: 53013
cat('# studies: ', trends[, length(unique(STUDY_ID))], '\n')
# studies: 332
cat('Data points: ', trends[, sum(nyrBT)], '\n')
Data points: 293973
trends[, table(REALM)]
REALM
Freshwater Marine Terrestrial
1025 48647 3341
trends[, table(taxa_mod)]
taxa_mod
All Amphibians Benthos Birds
1705 379 4679 13741
Fish Invertebrates Mammals Marine invertebrates/plants
28473 2996 525 206
Plant Reptiles
305 4
trends[, table(taxa_mod, REALM)]
REALM
taxa_mod Freshwater Marine Terrestrial
All 0 1702 3
Amphibians 2 0 377
Benthos 0 4679 0
Birds 0 11099 2642
Fish 1006 27467 0
Invertebrates 15 2901 80
Mammals 0 478 47
Marine invertebrates/plants 0 206 0
Plant 1 115 189
Reptiles 1 0 3
With all covariates
# the cases we can compare
apply(trends[, .(Jtutrend, REALM, tempave.sc, tempave_metab.sc, seas.sc, microclim.sc, temptrend.sc, mass.sc, speed.sc, lifespan.sc, consumerfrac.sc, endothermfrac.sc, nspp.sc, thermal_bias.sc, npp.sc, human.sc)], MARGIN = 2, FUN = function(x) sum(!is.na(x)))
Jtutrend REALM tempave.sc tempave_metab.sc seas.sc microclim.sc
53013 53013 49916 49916 49916 51834
temptrend.sc mass.sc speed.sc lifespan.sc consumerfrac.sc endothermfrac.sc
49916 52820 52689 51540 47534 53013
nspp.sc thermal_bias.sc npp.sc human.sc
53013 49371 52863 53013
i <- trends[, complete.cases(Jtutrend, temptrend.sc, tempave_metab.sc, REALM, seas.sc, microclim.sc, npp.sc, mass.sc, speed.sc, lifespan.sc, consumerfrac.sc, thermal_bias.sc)]
cat('Overall # time-series: ', sum(i), '\n')
Overall # time-series: 43585
cat('# studies: ', trends[i, length(unique(STUDY_ID))], '\n')
# studies: 250
cat('Data points: ', trends[i, sum(nyrBT)], '\n')
Data points: 222824
trends[i, table(REALM)]
REALM
Freshwater Marine Terrestrial
1008 39735 2842
trends[i, table(taxa_mod)]
taxa_mod
All Amphibians Benthos Birds Fish Invertebrates Mammals Plant
521 12 590 11803 27372 2567 518 200
Reptiles
2
trends[i, table(taxa_mod, REALM)]
REALM
taxa_mod Freshwater Marine Terrestrial
All 0 520 1
Amphibians 2 0 10
Benthos 0 590 0
Birds 0 9221 2582
Fish 993 26379 0
Invertebrates 12 2484 71
Mammals 0 477 41
Plant 1 64 135
Reptiles 0 0 2
Try combinations of
And choose the one with lowest AIC (not run: takes a long time)
# fit models for variance structure
fixed <- formula(Jtutrend ~ REALM + tempave_metab.sc + seas.sc + microclim.sc + npp.sc + temptrend_abs.sc +
mass.sc + speed.sc + lifespan.sc + consumerfrac.sc + thermal_bias.sc)
i <- trends[, complete.cases(Jtutrend, REALM, tempave_metab.sc, seas.sc, microclim.sc, npp.sc, temptrend_abs.sc,
mass.sc, speed.sc, lifespan.sc, consumerfrac.sc, thermal_bias.sc)]
mods <- vector('list', 0)
mods[[1]] <- gls(fixed, data = trends[i,])
mods[[2]] <- gls(fixed, data = trends[i,], weights = varPower(-0.5, ~nyrBT))
mods[[3]] <- gls(fixed, data = trends[i,], weights = varPower(0.5, ~ abs(temptrend)))
mods[[4]] <- lme(fixed, data = trends[i,], random = ~1|taxa_mod2, control = lmeControl(opt = "optim"))
mods[[5]] <- lme(fixed, data = trends[i,], random = ~1|STUDY_ID, control = lmeControl(opt = "optim"))
mods[[6]] <- lme(fixed, data = trends[i,], random = ~1|taxa_mod2/STUDY_ID, control = lmeControl(opt = "optim"))
mods[[7]] <- lme(fixed, data = trends[i,], random = ~1|STUDY_ID/rarefyID, control = lmeControl(opt = "optim"))
mods[[8]] <- lme(fixed, data = trends[i,], random = ~1|taxa_mod2/STUDY_ID/rarefyID, control = lmeControl(opt = "optim"))
mods[[9]] <- lme(fixed, data = trends[i,], random = ~temptrend_abs.sc | taxa_mod)
mods[[10]] <- lme(fixed, data = trends[i,], random = ~temptrend_abs.sc | STUDY_ID)
mods[[11]] <- lme(fixed, data = trends[i,], random = ~temptrend_abs.sc | taxa_mod2/STUDY_ID, control = lmeControl(opt = "optim"))
mods[[12]] <- lme(fixed, data = trends[i,], random = list(STUDY_ID = ~ temptrend_abs.sc, rarefyID = ~1)) # includes overdispersion. new formula so that random slope is only for study level (not enough data to extend to rarefyID).
mods[[13]] <- lme(fixed, data = trends[i,], random = list(taxa_mod2 = ~ temptrend_abs.sc, STUDY_ID = ~ temptrend_abs.sc, rarefyID = ~1)) # 30+ min to fit
mods[[14]] <- lme(fixed, data = trends[i,], random = ~1|STUDY_ID, weights = varPower(-0.5, ~nyrBT))
mods[[15]] <- lme(fixed, data = trends[i,], random = ~1|taxa_mod2, weights = varPower(-0.5, ~nyrBT))
mods[[16]] <- lme(fixed, data = trends[i,], random = ~1|taxa_mod2/STUDY_ID, weights = varPower(-0.5, ~nyrBT))
mods[[17]] <- lme(fixed, data = trends[i,], random = ~1|STUDY_ID/rarefyID, weights = varPower(-0.5, ~nyrBT))
mods[[18]] <- lme(fixed, data = trends[i,], random = ~1|taxa_mod2/STUDY_ID/rarefyID, weights = varPower(-0.5, ~nyrBT))
mods[[19]] <- lme(fixed, data = trends[i,], random = ~temptrend_abs.sc|STUDY_ID, weights = varPower(-0.5, ~nyrBT))
mods[[20]] <- lme(fixed, data = trends[i,], random = list(STUDY_ID = ~ temptrend_abs.sc, rarefyID = ~1), weights = varPower(-0.5, ~nyrBT))
mods[[21]] <- lme(fixed, data = trends[i,], random = list(taxa_mod2 = ~ temptrend_abs.sc, STUDY_ID = ~ 1), weights = varPower(-0.5, ~nyrBT))
mods[[22]] <- lme(fixed, data = trends[i,], random = list(taxa_mod2 = ~ temptrend_abs.sc, STUDY_ID = ~ 1, rarefyID = ~1), weights = varPower(-0.5, ~nyrBT))
mods[[23]] <- lme(fixed, data = trends[i,], random = list(taxa_mod2 = ~ temptrend_abs.sc, STUDY_ID = ~ temptrend_abs.sc), weights = varPower(-0.5, ~nyrBT)) # singular precision warning with lmeControl(opt = 'optim') and convergence error without
mods[[24]] <- lme(fixed, data = trends[i,], random = list(taxa_mod2 = ~ temptrend_abs.sc, STUDY_ID = ~ temptrend_abs.sc, rarefyID = ~1), weights = varPower(-0.5, ~nyrBT)) # singular precision warning with lmeControl(opt = 'optim') and convergence error without
mods[[25]] <- lme(fixed, data = trends[i,], random = ~1|taxa_mod2, weights = varPower(-0.5, ~abs(temptrend)))
mods[[26]] <- lme(fixed, data = trends[i,], random = ~1|STUDY_ID, weights = varPower(-0.5, ~abs(temptrend)))
mods[[27]] <- lme(fixed, data = trends[i,], random = ~1|STUDY_ID/rarefyID, weights = varPower(-0.5, ~abs(temptrend)))
mods[[28]] <- lme(fixed, data = trends[i,], random = ~1|taxa_mod2/STUDY_ID/rarefyID, weights = varPower(-0.5, ~abs(temptrend)))
mods[[29]] <- lme(fixed, data = trends[i,], random = ~temptrend_abs.sc|STUDY_ID, weights = varPower(-0.5, ~abs(temptrend)))
mods[[30]] <- lme(fixed, data = trends[i,], random = ~temptrend_abs.sc|taxa_mod2/STUDY_ID, weights = varPower(-0.5, ~abs(temptrend)), control = lmeControl(opt = "optim"))
mods[[31]] <- lme(fixed, data = trends[i,], random = list(STUDY_ID = ~ temptrend_abs.sc, rarefyID = ~1), weights = varPower(-0.5, ~abs(temptrend)))
mods[[32]] <- lme(fixed, data = trends[i,], random = list(taxa_mod2 = ~ temptrend_abs.sc, STUDY_ID = ~ temptrend_abs.sc, rarefyID = ~1), weights = varPower(-0.5, ~abs(temptrend)), control = lmeControl(opt = "optim")) # singular precision warning
aics <- sapply(mods, AIC)
minaics <- aics - min(aics)
minaics
which.min(aics)
Chooses the random slopes (temptrend_abs) & intercepts for STUDY_ID, overdispersion, and variance scaled to number of years. We haven’t dealt with potential testing on the boundary issues here yet.
world <- map_data('world')
ggplot(world, aes(x = long, y = lat, group = group)) +
geom_polygon(fill = 'lightgray', color = 'white') +
geom_point(data = trends, aes(rarefyID_x, rarefyID_y, group = REALM, color = REALM), size = 0.5, alpha = 0.4) +
scale_color_brewer(palette="Set1", name = 'Realm') +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
panel.background = element_blank(), axis.line = element_line(colour = "black"),
legend.key=element_blank(),
axis.text=element_text(size=16),
axis.title=element_text(size=20)) +
labs(x = 'Longitude (°)', y = 'Latitude (°)')
Mostly northern hemisphere, but spread all over. No so much in Africa or much of Asia.
Lines are ggplot smoother fits by realm.
Strong trends with temperature change, but trends are pretty symmetric around no trend in temperature, which implies warming or cooling drives similar degree of community turnover. Some indication of less turnover for larger organisms (mass) Higher turnover on land with higher seasonality? More turnover for shorter-lived organisms? No really clear differences among realms.
Violin plots
trends[temptrend <= -0.5, temptrendtext := 'Cooling']
trends[abs(temptrend) <= 0.1, temptrendtext := 'Stable']
trends[temptrend >= 0.5, temptrendtext := 'Warming']
trends[abs(rarefyID_y) < 35, latzone := 'Subtropics']
trends[abs(rarefyID_y) >= 35 & abs(rarefyID_x) < 66.56339, latzone := 'Temperate']
trends[abs(rarefyID_y) >= 66.56339, latzone := 'Polar']
trends[, latzone := factor(latzone, levels = c('Subtropics', 'Temperate', 'Polar'))]
p1 <- ggplot(trends, aes(temptrend, Jtutrend, color = REALM, fill = REALM, size = nyrBT)) +
geom_point(na.rm = TRUE, shape = 16, alpha = 0.1) +
geom_smooth(data=subset(trends, abs(temptrend) < 0.75), method = 'gam', formula = y ~ s(x, bs = "cs"),
na.rm = TRUE) +
scale_color_brewer(palette="Set1", name = 'Realm') +
scale_fill_brewer(palette="Set1", name = 'Realm') +
labs(x = 'Temperature trend (°C/yr)', y = 'Jaccard turnover') +
scale_size_continuous(range = c(1, 8), breaks = c(2, 5, 20)) +
guides(size = guide_legend(title = 'Years',
override.aes = list(linetype=0, fill = NA, alpha = 1))) +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
panel.background = element_blank(), axis.line = element_line(colour = "black"),
legend.key=element_blank(),
legend.key.size = unit(0.5,"line"),
axis.text=element_text(size=8),
axis.title=element_text(size=10))
p2 <- ggplot(trends[!is.na(temptrendtext), ], aes(temptrendtext, Jtutrend)) +
geom_violin(draw_quantiles = c(0.25, 0.5, 0.75), fill = 'grey') +
labs(x = '', y = 'Jaccard turnover') +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
panel.background = element_blank(), axis.line = element_line(colour = "black"),
legend.key=element_blank(),
axis.text=element_text(size=8),
axis.title=element_text(size=10))
p3 <- ggplot(trends[abs(temptrend) >= 0.5 & !is.na(latzone), ], aes(latzone, Jtutrend)) +
geom_violin(draw_quantiles = c(0.25, 0.5, 0.75), fill = 'grey') +
labs(x = '', y = '') +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
panel.background = element_blank(), axis.line = element_line(colour = "black"),
legend.key=element_blank(),
axis.text=element_text(size=7),
axis.title=element_text(size=10))
grid.arrange(p1, p2, p3, ncol = 2, layout_matrix = rbind(c(1,1), c(2,3)),
heights=c(unit(0.66, "npc"), unit(0.34, "npc")))
Average rates of turnover
trends[abs(temptrend) >= 0.5, .(mean(Jtutrend), sd(Jtutrend)/sqrt(.N))] # turnover per year for locations changing temperature
trends[abs(temptrend) < 0.1, .(mean(Jtutrend), sd(Jtutrend)/sqrt(.N))] # not changing temperature
trends[temptrend >= 0.5, .(mean(Jtutrend), sd(Jtutrend)/sqrt(.N))] # warming
trends[temptrend <= -0.5, .(mean(Jtutrend), sd(Jtutrend)/sqrt(.N))] # cooling
trends[abs(temptrend) >= 0.5 & abs(rarefyID_y) < 35, .(mean(Jtutrend), sd(Jtutrend)/sqrt(.N))] # tropics and sub-tropics
trends[abs(temptrend) >= 0.5 & abs(rarefyID_y) >= 35 & abs(rarefyID_y) < 66.56339, .(mean(Jtutrend), sd(Jtutrend)/sqrt(.N))] # temperate
trends[abs(temptrend) >= 0.5 & abs(rarefyID_y) >= 66.56339, .(mean(Jtutrend), sd(Jtutrend)/sqrt(.N))] # arctic
Hexagonal bins
require(ggplot2)
ggplot(trends[REALM == 'Terrestrial',], aes(temptrend, Jtutrend)) +
geom_hex() +
scale_fill_gradient(trans = "log")
ggplot(trends[REALM == 'Marine',], aes(temptrend, Jtutrend)) +
geom_hex() +
scale_fill_gradient(trans = "log")
ggplot(trends[REALM == 'Freshwater',], aes(temptrend, Jtutrend)) +
geom_hex() +
scale_fill_gradient(trans = "log")
i <- trends[, !duplicated(rarefyID)]; sum(i)
par(mfrow=c(5,3))
beanplot(rarefyID_y ~ REALM, data = trends[i,], what = c(1,1,1,1), col = c("#CAB2D6", "#33A02C", "#B2DF8A"), border = "#CAB2D6", ylab = 'Latitude (degN)', ll = 0.05)
beanplot(tempave ~ REALM, data = trends[i,], what = c(1,1,1,1), col = c("#CAB2D6", "#33A02C", "#B2DF8A"), border = "#CAB2D6", ylab = 'Temperature (degC)', ll = 0.05)
beanplot(tempave_metab ~ REALM, data = trends[i,], what = c(1,1,1,1), col = c("#CAB2D6", "#33A02C", "#B2DF8A"), border = "#CAB2D6", ylab = 'Metabolic Temperature (degC)', ll = 0.05, bw = 'nrd0') # nrd0 bandwidth to calculation gap
beanplot(seas ~ REALM, data = trends[i,], what = c(1,1,1,1), col = c("#CAB2D6", "#33A02C", "#B2DF8A"), border = "#CAB2D6", ylab = 'Seasonality (degC)', ll = 0.05)
beanplot(microclim ~ REALM, data = trends[i,], what = c(1,1,1,1), col = c("#CAB2D6", "#33A02C", "#B2DF8A"), border = "#CAB2D6", ylab = 'Microclimates (degC)', ll = 0.05)
beanplot(temptrend ~ REALM, data = trends[i,], what = c(1,1,1,1), col = c("#CAB2D6", "#33A02C", "#B2DF8A"), border = "#CAB2D6", ylab = 'Temperature trend (degC/yr)', ll = 0.05)
beanplot(mass_mean_weight ~ REALM, data = trends[i,], what = c(1,1,1,1), col = c("#CAB2D6", "#33A02C", "#B2DF8A"), border = "#CAB2D6", ylab = 'Mass (g)', ll = 0.05, log = 'y')
beanplot(speed_mean_weight +1 ~ REALM, data = trends[i,], what = c(1,1,1,1), col = c("#CAB2D6", "#33A02C", "#B2DF8A"), border = "#CAB2D6", ylab = 'Speed (km/hr)', ll = 0.05, log = 'y')
beanplot(lifespan_mean_weight ~ REALM, data = trends[i,], what = c(1,1,1,1), col = c("#CAB2D6", "#33A02C", "#B2DF8A"), border = "#CAB2D6", ylab = 'Lifespan (yr)', ll = 0.05, log = 'y')
#beanplot(consfrac ~ REALM, data = trends[i,], what = c(1,1,1,1), col = c("#CAB2D6", "#33A02C", "#B2DF8A"), border = "#CAB2D6", ylab = 'Consumers (fraction)', ll = 0.05, log = '') # too sparse
#beanplot(endofrac ~ REALM, data = trends[i,], what = c(1,1,1,1), col = c("#CAB2D6", "#33A02C", "#B2DF8A"), border = "#CAB2D6", ylab = 'Endotherms (fraction)', ll = 0.05, log = '') # too sparse
beanplot(Nspp ~ REALM, data = trends[i,], what = c(1,1,1,1), col = c("#CAB2D6", "#33A02C", "#B2DF8A"), border = "#CAB2D6", ylab = 'Number of species', ll = 0.05, log = 'y')
beanplot(thermal_bias ~ REALM, data = trends[i & !is.na(thermal_bias),], what = c(1,1,1,1), col = c("#CAB2D6", "#33A02C", "#B2DF8A"), border = "#CAB2D6", ylab = 'Thermal bias (degC)', ll = 0.05)
beanplot(npp ~ REALM, data = trends[i,], what = c(1,1,1,1), col = c("#CAB2D6", "#33A02C", "#B2DF8A"), border = "#CAB2D6", ylab = 'NPP', ll = 0.05)
beanplot(human ~ REALM, data = trends[i,], what = c(1,1,1,1), col = c("#CAB2D6", "#33A02C", "#B2DF8A"), border = "#CAB2D6", ylab = 'Human impact score', ll = 0.05)
Marine are in generally warmer locations (seawater doesn’t freeze) Marine have much lower seasonality. Marine and freshwater have some very small masses (plankton), but much of dataset is similar to terrestrial. Marine has a lot of slow, crawling organisms, but land has plants. Land also has birds (fast).
Try static covariates plus interactions of abs temperature trend with each covariate:
Except for thermal bias: interact with temperature trend (not abs)
summary(modonlyTtrend)
Linear mixed-effects model fit by REML
Data: trends[i, ]
Random effects:
Formula: ~temptrend_abs.sc | STUDY_ID
Structure: General positive-definite, Log-Cholesky parametrization
StdDev Corr
(Intercept) 0.04460504 (Intr)
temptrend_abs.sc 0.01904043 0.156
Formula: ~1 | rarefyID %in% STUDY_ID
(Intercept) Residual
StdDev: 0.001713667 0.2946268
Variance function:
Structure: Power of variance covariate
Formula: ~nyrBT
Parameter estimates:
power
-1.23577
Fixed effects: Jtutrend ~ abs(temptrend) * REALM
Correlation:
(Intr) abs(t) REALMM REALMT a():REALMM
abs(temptrend) -0.175
REALMMarine -0.928 0.162
REALMTerrestrial -0.928 0.162 0.861
abs(temptrend):REALMMarine 0.173 -0.990 -0.164 -0.161
abs(temptrend):REALMTerrestrial 0.172 -0.984 -0.160 -0.166 0.974
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-8.46254351 -0.30730773 0.08443652 0.54660915 6.61661896
Number of Observations: 49916
Number of Groups:
STUDY_ID rarefyID %in% STUDY_ID
317 49916
summary(modonlyTtrendJbeta)
Linear mixed-effects model fit by REML
Data: trends[i2, ]
Random effects:
Formula: ~temptrend_abs.sc | STUDY_ID
Structure: General positive-definite, Log-Cholesky parametrization
StdDev Corr
(Intercept) 0.05659393 (Intr)
temptrend_abs.sc 0.02520050 0.021
Formula: ~1 | rarefyID %in% STUDY_ID
(Intercept) Residual
StdDev: 1.693577e-06 0.3590015
Variance function:
Structure: Power of variance covariate
Formula: ~nyrBT
Parameter estimates:
power
-1.413558
Fixed effects: Jbetatrend ~ abs(temptrend) * REALM
Correlation:
(Intr) abs(t) REALMM REALMT a():REALMM
abs(temptrend) -0.160
REALMMarine -0.928 0.148
REALMTerrestrial -0.931 0.149 0.864
abs(temptrend):REALMMarine 0.158 -0.991 -0.150 -0.147
abs(temptrend):REALMTerrestrial 0.157 -0.986 -0.146 -0.151 0.976
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-7.8997139 -0.1696297 0.2037637 0.6704994 6.6449583
Number of Observations: 49916
Number of Groups:
STUDY_ID rarefyID %in% STUDY_ID
317 49916
summary(modonlyTtrendHorn)
Linear mixed-effects model fit by REML
Data: trends[i3, ]
Random effects:
Formula: ~temptrend_abs.sc | STUDY_ID
Structure: General positive-definite, Log-Cholesky parametrization
StdDev Corr
(Intercept) 0.05792109 (Intr)
temptrend_abs.sc 0.02134495 0.134
Formula: ~1 | rarefyID %in% STUDY_ID
(Intercept) Residual
StdDev: 0.007146559 0.306766
Variance function:
Structure: Power of variance covariate
Formula: ~nyrBT
Parameter estimates:
power
-1.217277
Fixed effects: Horntrend ~ abs(temptrend) * REALM
Correlation:
(Intr) abs(t) REALMM REALMT a():REALMM
abs(temptrend) -0.154
REALMMarine -0.914 0.141
REALMTerrestrial -0.925 0.143 0.845
abs(temptrend):REALMMarine 0.153 -0.989 -0.142 -0.141
abs(temptrend):REALMTerrestrial 0.152 -0.984 -0.139 -0.146 0.973
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-6.16836348 -0.32501087 0.06399723 0.53830958 5.96234033
Number of Observations: 48800
Number of Groups:
STUDY_ID rarefyID %in% STUDY_ID
276 48800
if(file.exists('temp/modonlyTtrendsimpreml.rds')){
modonlyTtrendsimpreml <- readRDS('temp/modonlyTtrendsimpreml.rds')
} else {
require(MASS) # for stepAIC
modonlyTtrendml <- update(modonlyTtrend, method = 'ML')
modonlyTtrendsimp <- stepAIC(modonlyTtrendml, direction = 'backward')
modonlyTtrendsimpreml <- update(modonlyTtrendsimp, method = 'REML')
saveRDS(modonlyTtrendsimpreml, file = 'temp/modonlyTtrendsimpreml.rds')
}
if(file.exists('temp/modonlyTtrendJbetasimpreml.rds')){
modonlyTtrendJbetasimpreml <- readRDS('temp/modonlyTtrendJbetasimpreml.rds')
} else {
require(MASS) # for stepAIC
modonlyTtrendJbetaml <- update(modonlyTtrendJbeta, method = 'ML')
modonlyTtrendJbetasimp <- stepAIC(modonlyTtrendJbetaml, direction = 'backward')
modonlyTtrendJbetasimpreml <- update(modonlyTtrendJbetasimp, method = 'REML')
saveRDS(modonlyTtrendJbetasimpreml, file = 'temp/modonlyTtrendJbetasimpreml.rds')
}
Start: AIC=-165090
Jbetatrend ~ abs(temptrend) * REALM
Df AIC
<none> -165090
- abs(temptrend):REALM 2 -164091
if(file.exists('temp/modonlyTtrendHornsimpreml.rds')){
modonlyTtrendHornsimpreml <- readRDS('temp/modonlyTtrendHornsimpreml.rds')
} else {
require(MASS) # for stepAIC
modonlyTtrendHornml <- update(modonlyTtrendHorn, method = 'ML')
modonlyTtrendHornsimp <- stepAIC(modonlyTtrendHornml, direction = 'backward')
modonlyTtrendHornsimpreml <- update(modonlyTtrendHornsimp, method = 'REML')
saveRDS(modonlyTtrendHornsimpreml, file = 'temp/modonlyTtrendHornsimpreml.rds')
}
Start: AIC=-146055.2
Horntrend ~ abs(temptrend) * REALM
Df AIC
<none> -146055
- abs(temptrend):REALM 2 -145334
summary(modonlyTtrendsimpreml)
Linear mixed-effects model fit by REML
Data: trends[i, ]
Random effects:
Formula: ~temptrend_abs.sc | STUDY_ID
Structure: General positive-definite, Log-Cholesky parametrization
StdDev Corr
(Intercept) 0.04460504 (Intr)
temptrend_abs.sc 0.01904043 0.156
Formula: ~1 | rarefyID %in% STUDY_ID
(Intercept) Residual
StdDev: 0.001713667 0.2946268
Variance function:
Structure: Power of variance covariate
Formula: ~nyrBT
Parameter estimates:
power
-1.23577
Fixed effects: Jtutrend ~ abs(temptrend) * REALM
Correlation:
(Intr) abs(t) REALMM REALMT a():REALMM
abs(temptrend) -0.175
REALMMarine -0.928 0.162
REALMTerrestrial -0.928 0.162 0.861
abs(temptrend):REALMMarine 0.173 -0.990 -0.164 -0.161
abs(temptrend):REALMTerrestrial 0.172 -0.984 -0.160 -0.166 0.974
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-8.46254351 -0.30730773 0.08443652 0.54660915 6.61661896
Number of Observations: 49916
Number of Groups:
STUDY_ID rarefyID %in% STUDY_ID
317 49916
summary(modonlyTtrendJbetasimpreml)
Linear mixed-effects model fit by REML
Data: trends[i2, ]
Random effects:
Formula: ~temptrend_abs.sc | STUDY_ID
Structure: General positive-definite, Log-Cholesky parametrization
StdDev Corr
(Intercept) 0.05659393 (Intr)
temptrend_abs.sc 0.02520050 0.021
Formula: ~1 | rarefyID %in% STUDY_ID
(Intercept) Residual
StdDev: 1.693577e-06 0.3590015
Variance function:
Structure: Power of variance covariate
Formula: ~nyrBT
Parameter estimates:
power
-1.413558
Fixed effects: Jbetatrend ~ abs(temptrend) * REALM
Correlation:
(Intr) abs(t) REALMM REALMT a():REALMM
abs(temptrend) -0.160
REALMMarine -0.928 0.148
REALMTerrestrial -0.931 0.149 0.864
abs(temptrend):REALMMarine 0.158 -0.991 -0.150 -0.147
abs(temptrend):REALMTerrestrial 0.157 -0.986 -0.146 -0.151 0.976
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-7.8997139 -0.1696297 0.2037637 0.6704994 6.6449583
Number of Observations: 49916
Number of Groups:
STUDY_ID rarefyID %in% STUDY_ID
317 49916
summary(modonlyTtrendHornsimpreml)
Linear mixed-effects model fit by REML
Data: trends[i3, ]
Random effects:
Formula: ~temptrend_abs.sc | STUDY_ID
Structure: General positive-definite, Log-Cholesky parametrization
StdDev Corr
(Intercept) 0.05792109 (Intr)
temptrend_abs.sc 0.02134495 0.134
Formula: ~1 | rarefyID %in% STUDY_ID
(Intercept) Residual
StdDev: 0.007146559 0.306766
Variance function:
Structure: Power of variance covariate
Formula: ~nyrBT
Parameter estimates:
power
-1.217277
Fixed effects: Horntrend ~ abs(temptrend) * REALM
Correlation:
(Intr) abs(t) REALMM REALMT a():REALMM
abs(temptrend) -0.154
REALMMarine -0.914 0.141
REALMTerrestrial -0.925 0.143 0.845
abs(temptrend):REALMMarine 0.153 -0.989 -0.142 -0.141
abs(temptrend):REALMTerrestrial 0.152 -0.984 -0.139 -0.146 0.973
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-6.16836348 -0.32501087 0.06399723 0.53830958 5.96234033
Number of Observations: 48800
Number of Groups:
STUDY_ID rarefyID %in% STUDY_ID
276 48800
coefs <- summary(modonlyTtrend)$tTable
coefs2 <- summary(modonlyTtrendJbeta)$tTable
coefs3 <- summary(modonlyTtrendHorn)$tTable
par(las = 1, mai = c(0.8, 2, 0.1, 0.1))
rows1 <- which(grepl('temptrend', rownames(coefs)))
plot(0,0, col = 'white', xlim=c(-0.02, 0.85), ylim = c(1,length(rows1)+0.3),
yaxt='n', xlab = 'Turnover per |°C/yr|', ylab ='')
axis(2, at = length(rows1):1, labels = c('Freshwater', 'Marine', 'Terrestrial'), cex.axis = 0.7)
abline(v = 0, col = 'grey')
for(i in 1:length(rows1)){
x = coefs[rows1[i], 1]
x2 = coefs2[rows1[i], 1]
x3 = coefs3[rows1[i], 1]
if(i>1){ # if an interaction, add the intercept
x = x + coefs['abs(temptrend)', 1]
x2 = x2 + coefs2['abs(temptrend)', 1]
x3 = x3 + coefs3['abs(temptrend)', 1]
}
se = coefs[rows1[i], 2]
se2 = coefs2[rows1[i], 2]
se3 = coefs3[rows1[i], 2]
points(x, length(rows1) + 1 - i, pch = 16)
lines(x = c(x-se, x+se), y = c(length(rows1) + 1 - i, length(rows1) + 1 - i))
points(x2, length(rows1) + 1.1 - i, pch = 16, col = 'light grey')
lines(x = c(x2-se2, x2+se2), y = c(length(rows1) + 1.1 - i, length(rows1) + 1.1 - i), col = 'light grey')
points(x3, length(rows1) + 1.2 - i, pch = 16, col = 'dark grey')
lines(x = c(x3-se3, x3+se3), y = c(length(rows1) + 1.2 - i, length(rows1) + 1.2 - i), col = 'dark grey')
}
legend('bottomright', col = c('black', 'dark grey', 'light grey'), lwd = 1, pch = 16,
legend = c('Jaccard turnover', 'Jaccard total', 'Horn-Morisita'))
summary(modTfull1simpreml)
Linear mixed-effects model fit by REML
Data: trends[i, ]
Random effects:
Formula: ~temptrend_abs.sc | STUDY_ID
Structure: General positive-definite, Log-Cholesky parametrization
StdDev Corr
(Intercept) 0.05008648 (Intr)
temptrend_abs.sc 0.01884335 0.56
Formula: ~1 | rarefyID %in% STUDY_ID
(Intercept) Residual
StdDev: 0.001194902 0.3037795
Variance function:
Structure: Power of variance covariate
Formula: ~nyrBT
Parameter estimates:
power
-1.246946
Fixed effects: Jtutrend ~ temptrend_abs.sc + REALM + tempave.sc + tempave_metab.sc + seas.sc + microclim.sc + mass.sc + speed.sc + lifespan.sc + endothermfrac.sc + nspp.sc + temptrend.sc + thermal_bias.sc + npp.sc + human.sc + temptrend_abs.sc:REALM + temptrend_abs.sc:tempave_metab.sc + temptrend_abs.sc:seas.sc + temptrend_abs.sc:microclim.sc + temptrend_abs.sc:mass.sc + temptrend_abs.sc:lifespan.sc + temptrend_abs.sc:nspp.sc + temptrend.sc:thermal_bias.sc + temptrend_abs.sc:npp.sc + temptrend_abs.sc:human.sc + REALM:human.sc + temptrend_abs.sc:REALM:human.sc
Correlation:
(Intr) tmpt_. REALMMr REALMTr tmpv.s tmpv_. ses.sc mcrcl. mss.sc spd.sc
temptrend_abs.sc 0.360
REALMMarine -0.918 -0.333
REALMTerrestrial -0.916 -0.325 0.835
tempave.sc 0.004 0.005 -0.014 -0.005
tempave_metab.sc 0.052 0.023 -0.037 -0.023 -0.289
seas.sc -0.036 0.000 0.053 -0.006 0.268 -0.084
microclim.sc -0.005 0.007 0.012 -0.016 0.049 0.063 0.171
mass.sc 0.022 0.000 -0.009 -0.020 0.006 -0.551 -0.045 -0.015
speed.sc 0.013 0.005 -0.025 0.022 -0.003 0.154 -0.072 -0.034 -0.128
lifespan.sc 0.032 0.024 -0.015 -0.006 -0.047 0.750 0.095 0.046 -0.807 0.246
endothermfrac.sc 0.145 0.007 -0.055 -0.222 0.123 -0.212 0.070 0.017 0.037 0.023
nspp.sc -0.026 -0.009 -0.002 -0.014 -0.042 -0.013 -0.032 -0.053 -0.161 -0.043
temptrend.sc -0.006 -0.013 0.006 0.007 -0.050 -0.018 -0.041 -0.014 0.010 -0.009
thermal_bias.sc 0.035 0.007 -0.047 -0.016 0.760 -0.053 -0.156 -0.135 0.043 0.016
npp.sc 0.001 -0.002 -0.003 0.015 -0.322 0.237 -0.255 -0.221 -0.009 -0.032
human.sc -0.167 -0.020 0.157 0.152 0.018 -0.005 0.009 0.017 0.005 0.004
temptrend_abs.sc:REALMMarine -0.341 -0.955 0.364 0.308 0.000 -0.009 0.011 -0.002 0.004 -0.001
temptrend_abs.sc:REALMTerrestrial -0.319 -0.896 0.299 0.337 -0.029 -0.021 -0.027 -0.024 0.008 -0.001
temptrend_abs.sc:tempave_metab.sc 0.025 0.081 -0.008 -0.016 0.103 0.365 0.049 0.080 -0.302 0.047
temptrend_abs.sc:seas.sc 0.019 -0.014 -0.010 -0.039 0.126 0.026 0.222 0.029 -0.029 -0.017
temptrend_abs.sc:microclim.sc -0.005 -0.036 0.001 -0.001 0.017 0.100 0.060 0.312 -0.036 0.015
temptrend_abs.sc:mass.sc 0.002 -0.001 0.003 -0.004 -0.008 -0.273 -0.037 -0.025 0.457 -0.030
temptrend_abs.sc:lifespan.sc 0.009 0.050 0.001 0.003 -0.011 0.367 0.050 0.022 -0.385 0.064
temptrend_abs.sc:nspp.sc -0.008 0.021 -0.002 0.001 -0.028 0.014 0.001 0.013 -0.070 0.029
temptrend.sc:thermal_bias.sc 0.009 0.019 -0.008 -0.007 0.067 0.033 0.000 -0.032 -0.002 0.008
temptrend_abs.sc:npp.sc 0.002 -0.054 -0.007 0.012 -0.115 0.095 -0.161 -0.163 0.013 0.000
temptrend_abs.sc:human.sc -0.070 0.160 0.066 0.063 -0.004 -0.009 0.005 0.006 0.009 -0.007
REALMMarine:human.sc 0.167 0.021 -0.158 -0.152 -0.020 0.004 -0.017 -0.024 -0.005 -0.003
REALMTerrestrial:human.sc 0.166 0.020 -0.157 -0.152 -0.032 0.007 -0.022 -0.006 -0.002 -0.001
temptrend_abs.sc:REALMMarine:human.sc 0.070 -0.158 -0.066 -0.063 0.001 0.009 -0.009 -0.009 -0.009 0.007
temptrend_abs.sc:REALMTerrestrial:human.sc 0.069 -0.160 -0.066 -0.063 0.012 0.007 -0.003 -0.005 -0.008 0.006
lfspn. endth. nspp.s tmptr. thrm_. npp.sc hmn.sc tm_.:REALMM tm_.:REALMT
temptrend_abs.sc
REALMMarine
REALMTerrestrial
tempave.sc
tempave_metab.sc
seas.sc
microclim.sc
mass.sc
speed.sc
lifespan.sc
endothermfrac.sc -0.054
nspp.sc 0.107 0.051
temptrend.sc -0.007 -0.002 0.013
thermal_bias.sc -0.075 0.028 -0.100 0.003
npp.sc 0.053 -0.129 -0.161 -0.013 -0.139
human.sc 0.008 0.006 -0.015 0.008 0.015 -0.029
temptrend_abs.sc:REALMMarine -0.011 0.005 -0.003 0.008 -0.003 -0.006 0.020
temptrend_abs.sc:REALMTerrestrial -0.017 0.019 0.000 0.015 -0.023 0.020 0.019 0.855
temptrend_abs.sc:tempave_metab.sc 0.385 0.027 -0.011 -0.049 0.108 0.003 -0.001 -0.048 -0.098
temptrend_abs.sc:seas.sc 0.035 0.053 0.001 -0.011 0.087 -0.178 -0.003 0.046 -0.082
temptrend_abs.sc:microclim.sc 0.025 -0.043 0.003 0.000 -0.026 -0.142 0.002 0.044 0.002
temptrend_abs.sc:mass.sc -0.381 0.023 -0.060 0.013 0.017 0.008 0.010 0.011 0.017
temptrend_abs.sc:lifespan.sc 0.466 -0.031 0.053 -0.010 -0.027 0.011 -0.004 -0.032 -0.038
temptrend_abs.sc:nspp.sc 0.076 -0.029 0.303 0.015 -0.051 -0.017 -0.011 -0.053 -0.039
temptrend.sc:thermal_bias.sc 0.003 0.008 -0.010 -0.402 -0.015 -0.003 -0.006 -0.017 -0.020
temptrend_abs.sc:npp.sc 0.010 -0.062 -0.034 -0.016 -0.044 0.366 -0.012 0.040 0.053
temptrend_abs.sc:human.sc -0.008 0.000 -0.003 0.007 -0.008 -0.012 0.409 -0.151 -0.139
REALMMarine:human.sc -0.009 -0.006 0.013 -0.008 -0.014 0.018 -0.998 -0.021 -0.019
REALMTerrestrial:human.sc -0.013 -0.011 0.018 -0.007 -0.020 0.027 -0.998 -0.020 -0.019
temptrend_abs.sc:REALMMarine:human.sc 0.008 0.000 0.000 -0.006 0.006 0.008 -0.408 0.150 0.139
temptrend_abs.sc:REALMTerrestrial:human.sc 0.008 0.002 0.004 -0.007 0.011 0.013 -0.406 0.151 0.138
t_.:_. tmptrnd_bs.sc:s. tmptrnd_bs.sc:mc. tmptrnd_bs.sc:ms.
temptrend_abs.sc
REALMMarine
REALMTerrestrial
tempave.sc
tempave_metab.sc
seas.sc
microclim.sc
mass.sc
speed.sc
lifespan.sc
endothermfrac.sc
nspp.sc
temptrend.sc
thermal_bias.sc
npp.sc
human.sc
temptrend_abs.sc:REALMMarine
temptrend_abs.sc:REALMTerrestrial
temptrend_abs.sc:tempave_metab.sc
temptrend_abs.sc:seas.sc 0.170
temptrend_abs.sc:microclim.sc -0.012 0.058
temptrend_abs.sc:mass.sc -0.626 -0.016 0.009
temptrend_abs.sc:lifespan.sc 0.829 0.045 0.021 -0.813
temptrend_abs.sc:nspp.sc 0.080 -0.056 -0.174 -0.224
temptrend.sc:thermal_bias.sc 0.072 -0.030 -0.036 -0.004
temptrend_abs.sc:npp.sc 0.046 -0.169 -0.213 -0.032
temptrend_abs.sc:human.sc -0.030 -0.018 0.038 0.049
REALMMarine:human.sc 0.001 -0.001 -0.005 -0.009
REALMTerrestrial:human.sc -0.001 -0.002 0.001 -0.008
temptrend_abs.sc:REALMMarine:human.sc 0.030 0.009 -0.037 -0.048
temptrend_abs.sc:REALMTerrestrial:human.sc 0.023 0.014 -0.027 -0.045
tmptrnd_bs.sc:l. tmptrnd_bs.sc:ns. tm.:_. tmptrnd_bs.sc:np.
temptrend_abs.sc
REALMMarine
REALMTerrestrial
tempave.sc
tempave_metab.sc
seas.sc
microclim.sc
mass.sc
speed.sc
lifespan.sc
endothermfrac.sc
nspp.sc
temptrend.sc
thermal_bias.sc
npp.sc
human.sc
temptrend_abs.sc:REALMMarine
temptrend_abs.sc:REALMTerrestrial
temptrend_abs.sc:tempave_metab.sc
temptrend_abs.sc:seas.sc
temptrend_abs.sc:microclim.sc
temptrend_abs.sc:mass.sc
temptrend_abs.sc:lifespan.sc
temptrend_abs.sc:nspp.sc 0.231
temptrend.sc:thermal_bias.sc -0.004 0.002
temptrend_abs.sc:npp.sc 0.042 -0.188 0.001
temptrend_abs.sc:human.sc -0.027 -0.032 -0.004 -0.034
REALMMarine:human.sc 0.003 0.010 0.005 0.008
REALMTerrestrial:human.sc 0.002 0.013 0.004 0.012
temptrend_abs.sc:REALMMarine:human.sc 0.027 0.028 0.005 0.024
temptrend_abs.sc:REALMTerrestrial:human.sc 0.021 0.033 0.006 0.030
tmptrnd_bs.sc:h. REALMM: REALMT: t_.:REALMM:
temptrend_abs.sc
REALMMarine
REALMTerrestrial
tempave.sc
tempave_metab.sc
seas.sc
microclim.sc
mass.sc
speed.sc
lifespan.sc
endothermfrac.sc
nspp.sc
temptrend.sc
thermal_bias.sc
npp.sc
human.sc
temptrend_abs.sc:REALMMarine
temptrend_abs.sc:REALMTerrestrial
temptrend_abs.sc:tempave_metab.sc
temptrend_abs.sc:seas.sc
temptrend_abs.sc:microclim.sc
temptrend_abs.sc:mass.sc
temptrend_abs.sc:lifespan.sc
temptrend_abs.sc:nspp.sc
temptrend.sc:thermal_bias.sc
temptrend_abs.sc:npp.sc
temptrend_abs.sc:human.sc
REALMMarine:human.sc -0.408
REALMTerrestrial:human.sc -0.408 0.997
temptrend_abs.sc:REALMMarine:human.sc -0.998 0.409 0.407
temptrend_abs.sc:REALMTerrestrial:human.sc -0.993 0.405 0.404 0.991
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-7.59387150 -0.34529883 0.05138559 0.55509022 6.91268058
Number of Observations: 43585
Number of Groups:
STUDY_ID rarefyID %in% STUDY_ID
250 43585
coefs <- summary(modTfull1simpreml)$tTable
par(las = 1, mai = c(0.5, 3, 0.1, 0.1))
rows1 <- which(!grepl('Intercept', rownames(coefs)))
plot(0,0, col = 'white', xlim=c(-0.02, 0.08), ylim = c(1,length(rows1)), yaxt='n', xlab = '', ylab ='')
axis(2, at = length(rows1):1, labels = rownames(coefs)[rows1], cex.axis = 0.7)
abline(v = 0, col = 'grey')
for(i in 1:length(rows1)){
x = coefs[rows1[i], 1]
se = coefs[rows1[i], 2]
points(x, length(rows1) + 1 - i, pch = 16)
lines(x = c(x-se, x+se), y = c(length(rows1) + 1 - i, length(rows1) + 1 - i))
}
resids <- resid(modTfull1simpreml)
preds <- getData(modTfull1simpreml)
col = '#00000033'
cex = 0.5
par(mfrow = c(4,4))
boxplot(resids ~ preds$REALM, cex = cex, col = col)
plot(preds$temptrend_abs.sc, resids, cex = cex, col = col)
plot(preds$temptrend.sc, resids, cex = cex, col = col)
plot(preds$tempave.sc, resids, cex = cex, col = col)
plot(preds$tempave_metab.sc, resids, cex = cex, col = col)
plot(preds$seas.sc, resids, cex = cex, col = col)
plot(preds$microclim.sc, resids, cex = cex, col = col)
plot(preds$mass.sc, resids, cex = cex, col = col)
plot(preds$speed.sc, resids, cex = cex, col = col)
plot(preds$lifespan.sc, resids, cex = cex, col = col)
plot(preds$consumerfrac.sc, resids, cex = cex, col = col)
plot(preds$endothermfrac.sc, resids, cex = cex, col = col)
plot(preds$nspp.sc, resids, cex = cex, col = col)
plot(preds$thermal_bias.sc, resids, cex = cex, col = col)
plot(preds$npp.sc, resids, cex = cex, col = col)
plot(preds$human.sc, resids, cex = cex, col = col)
i2 <- trends[, complete.cases(Jbetatrend, REALM, tempave.sc, tempave_metab.sc, seas.sc, microclim.sc,
temptrend.sc, temptrend_abs.sc, mass.sc, speed.sc, lifespan.sc,
consumerfrac.sc, endothermfrac.sc, nspp.sc, thermal_bias.sc, npp.sc, human.sc)]
i3 <- trends[, complete.cases(Horntrend, REALM, tempave.sc, tempave_metab.sc, seas.sc, microclim.sc,
temptrend.sc, temptrend_abs.sc, mass.sc, speed.sc, lifespan.sc,
consumerfrac.sc, endothermfrac.sc, nspp.sc, thermal_bias.sc, npp.sc, human.sc)]
randef <- list(STUDY_ID = ~ temptrend_abs.sc, rarefyID = ~1)
varef <- varPower(-0.5, ~nyrBT)
# full models
if(file.exists('temp/modTfullJbeta.rds')){
modTfullJbeta <- readRDS('temp/modTfullJbeta.rds')
} else {
modTfullJbeta <- lme(Jbetatrend ~ temptrend_abs.sc*REALM +
temptrend_abs.sc*tempave.sc +
temptrend_abs.sc*tempave_metab.sc +
temptrend_abs.sc*seas.sc +
temptrend_abs.sc*microclim.sc +
temptrend_abs.sc*mass.sc +
temptrend_abs.sc*speed.sc +
temptrend_abs.sc*lifespan.sc +
temptrend_abs.sc*consumerfrac.sc +
temptrend_abs.sc*endothermfrac.sc +
temptrend_abs.sc*nspp.sc +
temptrend.sc*thermal_bias.sc +
temptrend_abs.sc*npp.sc +
temptrend_abs.sc*human.sc*REALM,
random = randef, weights = varef, data = trends[i2,], method = 'REML')
saveRDS(modTfullJbeta, file = 'temp/modTfullJbeta.rds')
}
if(file.exists('temp/modTfullHorn.rds')){
modTfullHorn <- readRDS('temp/modTfullHorn.rds')
} else {
modTfullHorn <- lme(Horntrend ~ temptrend_abs.sc*REALM +
temptrend_abs.sc*tempave.sc +
temptrend_abs.sc*tempave_metab.sc +
temptrend_abs.sc*seas.sc +
temptrend_abs.sc*microclim.sc +
temptrend_abs.sc*mass.sc +
temptrend_abs.sc*speed.sc +
temptrend_abs.sc*lifespan.sc +
temptrend_abs.sc*consumerfrac.sc +
temptrend_abs.sc*endothermfrac.sc +
temptrend_abs.sc*nspp.sc +
temptrend.sc*thermal_bias.sc +
temptrend_abs.sc*npp.sc +
temptrend_abs.sc*human.sc*REALM,
random = randef, weights = varef, data = trends[i3,], method = 'REML')
saveRDS(modTfullHorn, file = 'temp/modTfullHorn.rds')
}
summary(modTfullJbeta)
summary(modTfullHorn)
# simplify
if(file.exists('temp/modTfullJbetasimp.rds')){
modTfullJbetasimp <- readRDS('temp/modTfullJbetasimp.rds')
} else {
require(MASS) # for stepAIC
modTfullJbetaml <- update(modTfullJbeta, method = 'ML')
modTfullJbetasimp <- stepAIC(modTfullJbetaml, direction = 'backward')
saveRDS(modTfullJbetasimp, file = 'temp/modTfullJbetasimp.rds')
}
if(file.exists('temp/modTfullHornsimp.rds')){
modTfullHornsimp <- readRDS('temp/modTfullHornsimp.rds')
} else {
require(MASS) # for stepAIC
modTfullHornml <- update(modTfullHorn, method = 'ML')
modTfullHornsimp <- stepAIC(modTfullHornml, direction = 'backward')
saveRDS(modTfullHornsimp, file = 'temp/modTfullHornsimp.rds')
}
summary(modTfullJbetasimp)
summary(modTfullHornsimp)
if(file.exists('temp/modTfullJbetasimpreml.rds')){
modTfullJbetasimpreml <- readRDS('temp/modTfullJbetasimpreml.rds')
} else {
modTfullJbetasimpreml <- update(modTfullJbetasimp, method = 'REML')
saveRDS(modTfullJbetasimpreml, file = 'temp/modTfullJbetasimpreml.rds')
}
if(file.exists('temp/modTfullHornsimpreml.rds')){
modTfullHornsimpreml <- readRDS('temp/modTfullHornsimpreml.rds')
} else {
modTfullHornsimpreml <- update(modTfullHornsimp, method = 'REML')
saveRDS(modTfullHornsimpreml, file = 'temp/modTfullHornsimpreml.rds')
}
# plot coefs
coefs2 <- summary(modTfullJbetasimpreml)$tTable
coefs3 <- summary(modTfullHornsimpreml)$tTable
varstoplot <- unique(c(rownames(coefs2), rownames(coefs3)))
rows1 <- which(!grepl('Intercept|REALM', varstoplot) | grepl(':', varstoplot)) # vars to plot in first graph
rows1_2 <- which(rownames(coefs2) %in% varstoplot[rows1]) # rows in coefs2
rows1_3 <- which(rownames(coefs3) %in% varstoplot[rows1]) # rows in coefs3
xlims1 <- range(c(coefs2[rows1_2,1] - coefs2[rows1_2,2],
coefs2[rows1_2,1] + coefs2[rows1_2,2],
coefs3[rows1_3,1] - coefs3[rows1_3,2],
coefs3[rows1_3,1] + coefs3[rows1_3,2]))
rows2 <- which(grepl('REALM', varstoplot) & !grepl(':', varstoplot)) # vars to plot in 2nd graph
rows2_2 <- which(rownames(coefs2) %in% varstoplot[rows2]) # rows in coefs2
rows2_3 <- which(rownames(coefs3) %in% varstoplot[rows2]) # rows in coefs3
xlims2 <- range(c(coefs2[rows2_2,1] - coefs2[rows2_2,2],
coefs2[rows2_2,1] + coefs2[rows2_2,2],
coefs3[rows2_3,1] - coefs3[rows2_3,2],
coefs3[rows2_3,1] + coefs3[rows2_3,2]))
cols <- c('black', 'grey') # for Jbeta and Horn models, respectively
offs1 <- 0.1 # offset vertically for the two models
offs2 <- 0.01 # offset vertically for the two models (plot 2)
par(mfrow=c(1,2), las = 1, mai = c(0.5, 3, 0.1, 0.1))
plot(0,0, col = 'white', xlim=xlims1, ylim = c(1,length(rows1)), yaxt='n', xlab = '', ylab ='')
axis(2, at = length(rows1):1, labels = varstoplot[rows1], cex.axis = 0.7)
abline(v = 0, col = 'grey')
for(i in 1:length(rows1)){
if(varstoplot[rows1[i]] %in% rownames(coefs2)){
x = coefs2[rownames(coefs2) == varstoplot[rows1[i]], 1]
se = coefs2[rownames(coefs2) == varstoplot[rows1[i]], 2]
points(x, length(rows1) + 1 - i + offs1, pch = 16, col = cols[1])
lines(x = c(x-se, x+se), y = c(length(rows1) + 1 - i + offs1, length(rows1) + 1 - i + offs1), col = cols[1])
}
if(varstoplot[rows1[i]] %in% rownames(coefs3)){
x = coefs3[rownames(coefs3) == varstoplot[rows1[i]], 1]
se = coefs3[rownames(coefs3) == varstoplot[rows1[i]], 2]
points(x, length(rows1) + 1 - i - offs1, pch = 16, col = cols[2])
lines(x = c(x-se, x+se), y = c(length(rows1) + 1 - i - offs1, length(rows1) + 1 - i - offs1), col = cols[2])
}
}
plot(0,0, col = 'white', xlim=xlims2, ylim = c(0.9, length(rows2) + 0.1), yaxt='n', xlab = '', ylab ='')
axis(2, at = length(rows2):1, labels = varstoplot[rows2], cex.axis = 0.7)
abline(v = 0, col = 'grey')
for(i in 1:length(rows2)){
if(varstoplot[rows2[i]] %in% rownames(coefs2)){
x = coefs2[rownames(coefs2) == varstoplot[rows2[i]], 1]
se = coefs2[rownames(coefs2) == varstoplot[rows2[i]], 2]
points(x, length(rows2) + 1 - i + offs2, pch = 16, col = cols[1])
lines(x = c(x-se, x+se), y = c(length(rows2) + 1 - i + offs2, length(rows2) + 1 - i + offs2), col = cols[1])
}
if(varstoplot[rows2[i]] %in% rownames(coefs3)){
x = coefs3[rownames(coefs3) == varstoplot[rows2[i]], 1]
se = coefs3[rownames(coefs3) == varstoplot[rows2[i]], 2]
points(x, length(rows2) + 1 - i - offs2, pch = 16, col = cols[2])
lines(x = c(x-se, x+se), y = c(length(rows2) + 1 - i - offs2, length(rows2) + 1 - i - offs2), col = cols[2])
}
}
Black is for Jaccard total turnover (pres/abs), grey is for Morisita-Horn turnover (considers abundance)